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Abstract -We study the Poisson-Boltzmann equation in the context of dense charged fluids 
where steric effects become important. We generalise the lattice gas theory by introducing a 
Flory-Huggins entropy for ions of differing volumes and then compare the effective free energy 
density to other approximations, valid for more realistic equations of state, such as the Carnahan- 
Starling approximation and find strong differences in the shapes of the free energy functions. We 
solve the Carnahan-Starling model in the high density limit, and demonstrate a slow, power-law 
convergence at high potentials. We elucidate how equivalent convex free energy functions can 
be constructed that describe steric effects in a manner which is more convenient for numerical 
minimisation. 


Introduction. — In the theory of ionic liquids jT] 
steric effects are of particular importance since the packing 
of ions can be especially dense [2j. The most common and 
simplest analytic approach to these effects is via the lat¬ 
tice gas mean-field approximation 3 . This methodology 
can be furthermore extended to a general local thermody¬ 
namic approach for any model of inhomogeneous fluids [dj. 
In this way one can connect the equation of state for any 
reference uncharged fluid, not only a lattice gas, with a 
full description of the same fluid with charged particles on 
the mean-field electrostatics level, generalizing in this way 
the Poisson-Boltzmann theory with consistent inclusion of 
packing effects. This approach is particularly relevant for 
analysis of dense electric double layers as arise in the con¬ 
text of ionic liquids or dense Coulomb fluids in general 5|. 
We will use this general local thermodynamics approach in 
conjunction with two such model equations of state: the 
asymmetric lattice gas approximation and the asymmetric 
Carnahan-Starling approximation. The size asymmetry as 
well as the charge asymmetry, Fig. [l] that this approach 
allows us to analyse, are fundamentally important for un¬ 
derstanding the nature of the electrostatic double layers. 

In what follows we will first formulate the local ther¬ 
modynamics mean-field approach to Coulomb fluids and 
then apply it to the asymmetric lattice gas, derived within 
the Flory-Huggins lattice approximation, comparing its re¬ 


sults with the asymmetric Carnahan-Starling approxima¬ 
tion. As a sideline we also derive several useful general 
relations valid specifically for the asymmetric lattice gas 
approximation in the context of electrostatic double lay¬ 
ers. 

General formulation. — We proceed by studying the 
Legendre transform of the free energy density /(ci, C 2 ) of 
an isothermal ( T = const) binary mixture 

/(ci,c 2 ) - M 1 C 1 - M 2 C 2 , (1) 


where C\ >2 are the densities of the two components, and 
the chemical potentials /xi s 2 are defined as 


Mi,2 — 


df(q,c 2 ) 

3ci, 2 


( 2 ) 


By the well known thermodynamic relationships [4] the 
Legendre transform eq. 0 equals 


/(ci,c 2 ) - 


d/(ci,c 2 ) <9 /(ci,c 2 ) 


dc 1 


Cl 


dc 2 


C 2 = -p(ci,c 2 ), (3) 


where p(ci,c 2 ) is the thermodynamic pressure, or the 
equation of state. For the inhomogeneous case we now 
invoke the local thermodynamic approximation so that 
the inhomogeneity is described solely via the coordinate 
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dependence of the densities, but the form of the thermo¬ 
dynamic potential remains the same as in the bulk, 

T = d 3 r (/(ci, c 2 ) - MiG - H 2 C 2 ) = - d 3 r p(pi, p 2 )- 
Jv Jv 

( 4 ) 

In the case of charged particles one needs to consider also 
the electrostatic energy and its coupling to the density 
of the particles via the Poisson equation, on top of the 
reference free energy of uncharged particles. The corre¬ 
sponding thermodynamic potential of the charged binary 
mixture then assumes the form 


of the binary components from eq. § , make a substitution 
Mi,2 —> Mi,2 =F ez^ 2 ip 

and finally derive the Euler-Lagrange equation for the lo¬ 
cal electrostatic potential of the form 

+ (8 ) 

dip 

which generalizes a form derived within a symmetric lat¬ 
tice gas approximation j8]. Invoking furthermore the 
Gibbs-Duhem relation 


J‘[ci,c 2 ,D] = / d 3 r (/(ci, c 2 ) -MiG - M2C2) + 

Jv 

+ Jd 3 r — ip(V • D - e{z 1 c 1 - z 2 c 2 ))j , (5) 

where D = D(r) is the dielectric displacement field, 
£ = eeo with e the relative dielectric permittivity, Ziy are 
the valencies of the two charged species and ip = ip( r) is 
now the Lagrange multiplier field that ensures the local 
imposition of Gauss’ law 6|. We can write this expression 
in an alternative form as 

T = / d 3 r (/(ci, c 2 ) - (mi - ez\ip)ci - (M2 +ez 2 ip)c 2 ) + 

Jv 

+ If r {^-* VD )- (6) 

Invoking now the identity eq. @ , discarding the boundary 
terms and minimizing with respect to D, we get the final 
form of the inhomogeneous thermodynamic potential 


dp 

ci ,2 = 7 ;- 

C^Ml.2 

we derive the Poisson equation as 
dp{p 1 - eziip, p 2 + ez 2 ip) dp dp 

-a*- = -“ i s 

= -e(zici - z 2 c 2 ) = q. (9) 

where q is the local charge density. Note that the charge 
density is a derivative w.r.t. potential of a single func¬ 
tion, a simple test of consistency for any proposed theory. 
Together with eq. 0 this constitutes a generalisation of 
the Poisson-Boltzmann theory for any model of the fluid 
expressible via an equation of state in the local thermody¬ 
namic approximation. This also generalizes some results 
previously derived only for the lattice gas. 

In the case of a single or two planar surfaces, with a 
normal in the direction of the 2 -axis, so that ip( r) = ip(z), 
the Poisson-Boltzmann equation possesses a first integral 
of the form 


= -/ d 3 r (Jeiyip ) 2 + p(mi - eziip, M2 + ez 2 tp)) ■ 

Jv 

( 7 ) 

In the case of charged boundaries one needs to add a sur¬ 
face term <f> s ipD n dS , where D n is the normal component 
of the electric displacement field at the surface, to the 
above equation. While the derivation of eq. 0 proceeded 
entirely on the mean-field level, it can be extended to the 
case when the Coulomb interactions are included exactly 
and the mean potential becomes the fluctuating local po¬ 
tential in a functional integral representation of the parti¬ 
tion function [7|. 

Let us note that the signs of the electrostatic terms in 
eq. 0 are consistent with the definition of the grand 
canonical partition function, i.e. Q = — pV, with 

fi(A,0) = E ^=1 \ n Q{N,P)/N\, where Q(N,f3) is the 
canonical partition function for N particles. The abso¬ 
lute activity is defined as A = Since electrostatic 

interactions enter with a Boltzmann factor, A = —> 

g/3/j'Tez, , 21 /y where — is valid for positive and + for negative 
ions. 

For any equation of state p(p 1 , p 2 ) or indeed any model 
free energy f(ci,c 2 ) of the reference uncharged system, 
one now needs to evaluate the proper chemical potentials 


| eip' 2 (z) - p(m 1 - ez y ip(z), p 2 + ez 2 ip{z)) = -p 0 (10) 

where po is an integration constant equal to the osmotic 
pressure of the ions and determined by the boundary 
conditions. The disjoining (interaction) pressure for two 
charged surfaces, II, is then obtained by subtracting the 
bulk contribution from the osmotic pressure po- The first 
integral of the Euler-Lagrange equation can be used to 
construct an explicit ID solution, ip = ip(z) by quadra¬ 
ture. 

In the limiting case of an ideal gas, with the van’t Hoff 
equation of state p(ci,c 2 ) = fc_eT(ci + C 2 ) it is straight¬ 
forward to see that the above theory reduces exactly to 
the Poisson-Boltzmann approximation |9j. Furthermore, 
for the binary, symmetric lattice-gas 

k T 

p{ci, c 2 ) = —log (1 — a 3 (ci + c 2 )) (11) 

where a is the cell size 1(5], the above formalism yields 
the results discussed at length by Kornyshev [3 . From 
eq. we also see one of the weaknesses of the lattice 
gas approach, as the pressure diverges only very weakly 
at close packing. We will compare with a more realistic 
equation of state later in this paper. 
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Figure 1: A schematic representation of the asymmet¬ 
ric system composed of particles with unequal sizes and 
unequal valences in the vicinity of a charged surface. 


Asymmetric lattice gas. — We start with the free 
energy density of mixing for a three component lattice 
gas system composed of species ”1” at concentration ci, 
itself composed of N\ subunits, and species ”2” at con¬ 
centration C 2 , itself composed of A 2 subunits, in a solvent 
of (water) molecules of diameter a. It can be expressed 
rather straightforwardly in terms of the volume fractions 
0i, 0 2 after realizing that it is equivalent to the problem of 
polydisperse polymer mixtures on the Flory-Huggins lat¬ 
tice level 11 . For a two component system the free energy 


of mixing can be derived simply as 12 


/(</> i,h) a 3 </>i . , , 02 . , . 

--= AT lo § 0 1 + AT lo S 0 2+ 

k B T N x N 2 

+ (1 - 0i. - 02 ) log (1 - 4>i - 0 2 ), 
where the volume fractions <j >\, 0 2 are defined as 
01,2 = = R 3 2 c 1,2i 


( 12 ) 


(13) 


and Nip = (^ 1 , 2 /a) measures the relative volumes of 
species 1 and 2, with radii 7?i, 2 , compared to the solvent 
with radius a. While the size-symmetric lattice gas has a 
venerable history (for an excellent review see Ref. [1]) there 
have been fewer previous attempts to master the lattice 
gas mixtures in the context of size-asymmetric electrolytes 
[5 13-16 and the simple connection with the entropy of 


lattice polymers has apparently not been noted before. 
The chemical potential is then obtained as 


Mi ,2 = 


df(d,c 2 ) 5/(</>i, <j> 2 ) 


dci , 2 


50 


a 3 A,; 


1,2 


(14) 


that can be evaluated explicitly yielding 


Ml ,2 = log 01,2 + 1 - N 1>2 (log (1 - 01 - 02) + 1) ■ (15) 

The Legendre transform eq. (jT]) then yields the osmotic 
pressure, again as a function of both volume fractions 


p(0i,0 2 ) a 3 


k B T 

1 - 


= log (1 - 01 - 02 )T 


A, 


1 - 


1 

A0 


(16) 


The form of this result is revealing as it states that the 
osmotic pressure is basically the lattice gas pressure of a 
symmetric mixture, corrected by the fact that Ai, 2 sub¬ 
units of the species ”1” and ”2” do not represent separate 
degrees of freedom. Obviously, for a symmetric system 
with Ai, 2 = 1 this reduces exactly to the lattice gas sym¬ 
metric binary mixture expression, eq. 


Introducing = Mi, 2 +-^ 1,2 _ 1 we can rewrite eq. (15) 


as 


* 1,2 = (1 - 01 - 0 2 ) JVl ' 2 e Ml ’ 2 


(17) 


Using this relation we can derive an explicit equation for 


= (1 - 01 - 02) 


(18) 


of the form 


u(l + u Nl ~ 1 e iil + u^-V 2 ) = 1, (19) 

that yields u = u(p,i, p, 2 ; A/, N 2 ). This allows us to fi¬ 
nally write the osmotic pressure as a function of the two 
densities 

- P ^ Cl ' Cl r l a = log (1 - a 3 (c 1 N l + c 2 N 2 ))+ 
k B T 

+ a 3 ci (Ai - 1) + a 3 c 2 (N 2 - 1), (20) 


or of the two chemical potentials through u = 
m(mi,M 2 ! A 1; N 2 ) as 


p(Mi,M2) a 3 


k B T 


= log u+ 


u Nl e' 11 


1 - 


1 

N, 


■ u N2 e^ 2 


1 - 


A 2 


( 21 ) 


In the case of ions of the same size, we can set without 
any loss of generality, that Ai = A 2 = 1, so that 


p(Mi,M 2 ) « 3 

k B T 


= log u{jh,jl 2 ) = — log (l + e Ml + e^ 2 ), 

( 22 ) 


a standard expression for the symmetric lattice gas [l7] . 

Above equations present a complete set of relations sat¬ 
isfied by the asymmetric lattice gas, being a mixture of two 
differently sized ions. The addition of mean-field electro¬ 
static interactions eq ([8]) then modifies solely the chemical 
potentials so that 


p(M i,M 2 )—t p(mi — eziip, M 2 + ez 2 tp) (23) 


if the two species are oppositely charged, which we as¬ 
sume. The corresponding Poisson-Boltzmann equation is 
then obtained from eq. [8] and eq. [9] in the form 


eV 2 0 = —e (zi5 Ml - z 2 d^ 2 )p(pi - ez^tp, p 2 + ez 2 ip) = 

= ( 24 ) 


where 0 i, 2 ( 0 ) are obtained from eq. 
Mi — eziip, fi 2 —> jl 2 + ez 2 ip. 
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In complete analogy with the case of polyelectrolytes 
with added salt [18] it is clear that electroneutrality of the 
asymmetric lattice gas in the bulk is achieved only if it is 
held at a non-zero electrostatic potential, tpo, that can be 
obtained from eq. (24) in an implicit form 


(iVi - N 2 )logu(tp 0 ) 


-(A i 


fa) + log 


N 1 Z 2 
N 2 zi ' 


(25) 


In what follows we then simply displace the origin of the 
electrostatic potential by ipo, the Donnan potential, inter¬ 
preted as the change in electrostatic potential across the 
bulk reservoir - ionic liquid interface, or equivalently as a 
Lagrange multiplier for the constraint of global electroneu¬ 
trality 


19 


Asymptotic behaviour of the lattice gas model. 

— We now consider the forms of the general equations 
derived above in the limiting cases of small and large elec¬ 
trostatic potential of the lattice gas model: 

Small potential and screening length. In the limit of 
ip —¥ 0 , one can derive 


p(pi + ezxip , fM 2 - ez 2 ip) = p(p!,n 2 ) + 

-e(zid Ml - z 2 d fl2 )p(p lt ii 2 )ip + 

+ ±e 2 (zid^ - z 2 d^) 2 p(p 1 ,p 2 ) ip 2 + Opip 3 ) (26) 


where we took into account eq. 0- Just as in the full non¬ 
linear case, see above, the term linear in ip is connected 
with the displaced electrostatic potential, eq. (25). The 
linearized form of ipo is then obtained approximately as 


straightforwardly that 

2 e 2 ^ Cf + 2 -i-2 

k 2 = - - - 

£ Q2f Q2f _ 

dc\ 802 


d 2 f , .2 a 2 / 

dcidci ' 1 dc '2 



. 4 n f„ N ,No u{z2lCl + Z2C2) + a3ciC2 ( ZlN2 + z ^ 2 

B 1 2 (1 + (Ni — l)a 3 NiCi + (N 2 — l)a 3 N 2 c 2 ) ’ 

(30) 


where we introduced the Bjerrum length l B = 
e 2 /(‘{-KeksT). In general the Debye length is therefore 
not a linear function of the concentrations. For the sym¬ 
metric lattice gas, Ni = N 2 = 1, and taking into ac¬ 
count the definition eq. [l8j the above result reduces 
to k 2 = 47 t£ b ( z i c i + z 2 c 2 ~ a 3 (zici — z 2 c 2 ) 2 ), which for 


bulk electroneutrality reduces further to the standard De¬ 


bye expression 17 


Large potential and close packing. The limits for ip —► 
±00 of a lattice gas can be derived as 


_ f e -(Ai+ezi</0/lVi ^ +QO 

U ~ \ e -(r2-ez 2 ^)/N 2 ^ _ 0Q 

implying 


(31) 


<Pl,2(lp ->• OO) 


1 

e -(p.i+ez 1 i/j)N2/N 1 +(p.2-ez2il’) 


(32) 


and 


<pi, 2 (ip -t - 00 ) 


e -(p,2-ez2it>)N 1 /N 2 + (p,i-ez 1 ip) 

1 


(33) 


ejzxd^ - z 2 d^ 2 )p(p 1 ,p 2 ) 
e 2 (z 1 d lll - z 2 d IJ2 ) 2 p{p 1 , p 2 ) 


(27) 


Obviously the expansion of the pressure for small values 
of the electrostatic potential is then quadratic in the dif¬ 
ference ip — Ipo- 

Furthermore, the Hessian of the pressure p(p 1 ,^ 2 ) is 
positive definite, i.e. 


where the upper formula is for ”1” and the lower for ”2”. 
Thus 


p(ip —> 00 ) 
k B T 
1 


= -(fa + ez 1 ip)/N 1 + 


_|_ [ l __ j _|_ e -{izi+ezixl>)N2/Nx e tp.2-ez2il>) j ^ 


1 

A 2 


(34) 


e 


2 


e 



2zi z 2 


d 2 p 
djM p 2 



= K 2 > 0 , 


(28) 


while from eq. 0 it follows that n is nothing but the in¬ 
verse Debye length expressed through the second deriva¬ 
tives of the pressure with respect to the chemical potentials 
of both charged species. Since the curvature tensor of the 
Legendre transform is the inverse of the curvature tensor 
of the function itself [20 , we can write 

V d 2p(M) d 2/(ci ’ cz) = Slk (29) 

^ dmOpm dCmdCk K ’ 


where all the matrices are 2x2. From here it follows rather 


and 


p(tp —> — 00 ) a 3 . ,, r 

-r 77 T-= -(^2 - ez 2 ip) N 2 + 

k B T 

1 p -(h2-ez2tl>)N 1 /N 2 Jpi+ez 1 ip) (i _;L \ i / 1_)L | 

V nJ v n 2 ) 

(35) 

The most striking feature of these limits is the linear 
behaviour of p(ip) for large positive or negative potentials, 
which gives rise to a V-like curve for symmetric particles. 
This linear behaviour is linked to the saturation of close 
packing of the lattice particles against a high potential 
surface. For particle of unequal size the two branches of 
p(ip) have different slopes. 
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A dense two component lattice gas. — For some 
cases in the theory of ionic liquids one can assume dense 
packing, without any intervening solvent, so that 0 i +02 = 
1. The corresponding free energy can then be cast into a 
simplified form 


/(0i,0 2 = 1 - 0i) N 2 a 3 
k B T 


0 i 

M 


log 01 + 02 log 02- 


(36) 


where M = N 1 /N 2 is the effective size of the species ”1” 
compared to species ”2”. This implies furthermore that 


df_ 

50i 


df_ 

<902 


= 0 


Pi T M pi 2 — 0. 


(37) 


The equation analogous to eq. (191 then assumes the form 

0r = (1 - 0 1 ) A V 1 " 1+M (38) 


and the Legendre transform of the free energy density fol¬ 
lows as 


p( 0 r ,02 = 1 - 0i) N 2 a 3 
k B T 


log 02 + 01 



(39) 


The Poisson-Boltzmann equation is then cast into a sim¬ 
plified form 

eV 2 0 = - e (zid^ - z 2 d t _ l2 )p{n 1 - eziip, p 2 + ez 2 i[>) = 

= AT (( Zl + Mz 2 )<t>i - Mz 2 ) , (40) 

iVi 

and the charge density is constrained to be between — 
and 

Asymmetric Carnahan-Starling approximation. 

— In order to show the interest and generality of our local 
thermodynamic approach we will now apply it in the case 
of the Carnahan - Sterling approximation for asymmetric 
binary hard sphere mixtures [21 . For a bulk, uncharged 
hard sphere fluid the Carnahan-Starling approximation is 
’’almost exact”. 

The excess pressure in the Carnahan-Starling approxi¬ 
mation derived via the ’’virial equation” is then equal to 




Figure 2: Top: the function g^(q) for a Carnahan-Starling 
fluid. It is an effective potential at imposed local charge 
density. The divergences correspond to the close packing 
of the fluid. Blue: symmetric particle volumes. Gold: 
asymmetric fluid with larger negative particles. This 
function is important in dual convex formulations of the 
Poisson-Boltzmann equation, see eq. (621. Bottom: the 
function p('ijj) for the same two sets of fluid parameters. 
Blue: Seen from afar the function displays a characteristic 
V-like behaviour for symmetric particles. The larger par¬ 
ticles, gold, give a smaller slope in the pressure function. 
The central part of the figure for symmetric particles is 
examined more closely in figure [3] 


with 


A \/£i£2 {Ri - R2) 2 

A ”- 1 — mr- 

The excess free energy then follows as 

feX ck B T 2) = ~i ( 1 -yi + y2 + y 3 ) + 

3y 2 + 2 y 3 3 1 - yi - U 2 - 

i -£ 5 ( i -£) 2 

+ ( y 3 - 1 ) log (1 £)? 


3 2/3 


(44) 


(45) 


Pexc(Cl,C 2 ) = (1 + £ + $ 2 ) - 3g(yi + &/ 2 ) - 3£ 3 ?/3 

ck B T (1 - £) 3 

where Ci ,2 are the densities of the two components and 


It is now straightforward to obtain the chemical potentials 
from the free energy as pi^ 2 = fJ,i^ 2 {ci, c 2 ), invert them and 
then obtain the pressure equation as p = p{pi,p, 2 ). 


4tt 3 

SI ,2 = —iti i 2 Ci , 2 


and 


? = £1 + £ 2 , 


(42) 


where Ri^ 2 are the hard sphere radii of the two species. 
Furthermore 


£i9?2 + £ 27 ?! 



Asymptotic behaviour for the Carnahan-Starling 
free energy density. — In the limit of large poten¬ 
tials 0, as occurs near an electrode, the second, wrongly 
charged, component of the fluid is excluded and the dom¬ 
inant physics is the packing of a single component system 
under the constraints coming from the electrostatic inter¬ 
actions. In this limit of £ —> 1 we can substitute y\^ 2 = 0 
and </3 = 1. The most important divergence in this limit 
(43) thus stems from the denominator in eq. (45). 

The free energy density of a Carnahan-Starling liquid 
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Figure 3: Comparison of symmetric Carnahan-Starling 
and lattice gas forms near ip = 0. Top p(ip), bottom charge 
density —q(ip) = dp/dip. Carnahan-Starling in blue, con¬ 
verges very slowly to the slope final value, (see figure. [4]), 
and is much broader at minimum compared to the lattice 
gas curve (red). Physical parameters are such that the 
blue and red curves have asymptotically the close packing 
charge density; the chemical potentials have been tuned to 
give the same Debye length in the bulk (corresponding to 
the same curvature at the origin in the curve for p). De¬ 
spite this double matching the curves are very different in 
functional form. The curve q(ip) for the Carnahan-Starling 
fluid converges extraordinarily slowly compared with the 
lattice gas approximation. 



Figure 4: Study of the evolution of — q = dp(ip)/dip in 
variables adapted to the Carnahan-Starling fluid, (blue). 
The large field behaviour is linear in 1/ip 1 ^ 3 , demonstrat¬ 
ing the correctness of arguments leading to eq. (48). Same 
data as fig. [3] for large positive potentials. The lattice gas 
model (red) gives much more rapid cross-over to satura¬ 
tion at large potentials. Horizontal line is to guide the eye 
and corresponds to the close packing charge density. 


near close packing has a singularity of the form 


/(c) 


c 0 k B T 
(1 - c/c 0 ) 2 


(46) 


Co is the close packing volume fraction of the component 
dominating near the electrode. With this assumption we 
can take the Legendre transform of the most singular, di¬ 
verging part of the free energy to find the large poten¬ 
tial limit of p(ip). For large positive ip (assuming that 
ez 2 ip p) this limit turns out to be 

p(ip) = z 2 ec 0 ip - ^(2c 0 k B T) 1/3 (z 2 ecoip) 2/3 (47) 

where we have used the fact that negative ions of valence 
z 2 dominate. 

In the high packing limit p{ip) is therefore linear in the 
potential, and is given by the spatial charge density at 
close packing exactly like for the lattice gas. However 
unlike the lattice gas the approach to the high field limit 
is very slow 

- q(iP) = z 2 ec 0 - 1 /3 (z 2 eco) 2/3 (48) 

The spatial charge density is negative for large positive 
potentials. The correctness of this law is demonstrated 
in figure [4] which plots (1/ coZ 2 e)dp/dip as a function of 
ip -1 / 3 . The curve linearly extrapolates to unity for large 
ip. There is a very clear contrast with the case of the lattice 
gas model where the cross-over to close packing occurs for 
much smaller values of the potential. 

Solution for the high field Carnahan-Starling limit. 
The solution for the generalized Poisson-Boltzmann equa¬ 
tion in the high field limit can be found from the solution 
of the integral problem 


ec 0 z 2 ip - p(ec 0 z 2 ip) 2/3 




(49) 

with /3 = 3(2cofcsT) 1 ^ 3 /2; we neglect p compared to 
z 2 eip. This integral can be transformed by substituting 
{ec 0 z 2 ip) 1/3 = y , giving 


dy y 2 

(■ V 3 - Py 2 + Po) 1/2 

a form which can be solved by using elliptic functions. If 
we make the further approximation that po is small we can 
find much simpler expressions: 

z{ip) - z 0 =^-^-[{(z 2 ecoip) 1/3 - P) / + 
ecoZ 2 

3P{{z 2 ec 0 ip) 1/3 - P) / } (51) 

Here, z(ip) gives the distance from a plate which corre¬ 
sponds to a potential ip. It is obviously the inverse function 



2 e 2 cgz| 


dz, 


(50) 
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Figure 5: Carnahan-Starling fluid. Blue: exact numerical 
calculation of the function z(tp). The full curve is given in 
the inset. The high field limit studied in the main panel. 
Red: approximation valid for large fields from eq. (51). 
Green profile calculated from assuming perfect packing of 
the fluid against electrode. Blue and red curves almost 
perfectly correspond, with no extra fitting parameters. 


of We perform a “numerically exact” calculation of 

the curve ip(z) in fig. ([5]) in the inset, where we place a 
positve electrode at z = 0. The main figure of fig. ^ 
contains three curves: The blue curve explodes part of the 
inset and is overlayed with a red curve corresponding to 
eq. (51). On this scale the results are indistinguishable. 
The green curve is evaluated by assuming perfect packing 
of the fluid against the electrode. Eq. (51) is clearly a 
much better description of the high electrostatic potential 
physics. 

Eq. (51) can also be combined with eq. (48) to find z(q) 
and thus the evolution of the spatial charge density with 
distance from an electrode as well as the variation of the 
local charge density with the potential. 


Differential Capacitance. — Together with the 
boundary condition D n = a, where cr is the surface charge 
density one can derive the equivalent of the Grahame equa¬ 
tion in the form 

-- p(Hi - eziipo, + ez2ipo) = -po, (52) 

Z£ 

assuming that the bounding surface is located at z = 0, 
i.e. x/jo = xp(z = 0). From the Grahame equation one can 
next derive the differential capacitance C as 


C(ipo) 


da(4>o) ee ( dp dp \ 

dil’o ct(V’o) V Zl d^ +Z2 dp 2 ) 

—ee(z 1 c 1 - z 2 c 2 ) 
±\/2e(p-po) 


(53) 


with ± depending on the sign of the surface charge. Taking 
into account the definition of the Bjerrum length, £b 

(£(*,) = 2 nk B Tt B (54) 

Wo) - Po 

Invoking the Poisson-Boltzmann equation for this case, an 
alternative form of the differential capacitance is 


C(V’o) =e(log^o)' = 


d 


=V2e-^—^p(pi - ezii.’ 0 ,p, 2 
94 ’ o 


ez 2 ^ 0 )-p 0 , (55) 


the form that we use in our numerical work. It is interest¬ 
ing to note that even if we shift the minimum of the curve 
p(i[b) to occur at 4 > = 0 this does not imply that ^ = 0 is 
also a stationary value of the differential capacitance. This 
is clearly visible in the curves of fig. © where in denser 
fluids the maximum of the curves is shifted to positive po¬ 
tentials. We mark the position of the minimum in p(4>) 
by a slight break in the solid lines. This displacement of 
the maximum of the capacitance from the minimum of p 
is trivially understood if one assumes that the expansion 
of p(4>) includes a term in -0 3 . We see that the qualitative 
behaviour of the curves generated for the lattice model, as 
well as the Carnahan-Starling fluid are rather similar. 


Conclusions. — By using general arguments based 
upon local thermodynamics, we generalized the Poisson- 
Boltzmann mean-field theory of Coulomb fluids to the case 
where the reference, uncharged fluid need not be ideal. We 
formulated the general theory in the particular cases of 
an asymmetric lattice gas and an asymmetric Carnahan- 
Starling liquid that describe steric effects at various levels 
of approximations and are particularly relevant for anal¬ 
ysis of dense electric double layers as arise in the context 
of ionic liquids or dense Coulomb fluids. Use of properties 
of Legendre transforms allows us to efficiently translate 
between forms of the free energy; this includes a standard 
formulation in terms of the electrostatic potential, and a 
dual formulation (see appendix) in terms of the electric 
displacement held. 

We analyzed in detail the size and charge asymmetry 
and their respective effects on the salient properties of 
electric double layers. As part of our analysis we also for¬ 
mulated an exact thermodynamic description of an asym¬ 
metric lattice gas, derived within the Flory-Huggins lattice 
approximation. This allows the lattice gas approximation, 
which in its symmetric form already serves as the most 
popular description of the steric effects in the context of 
the Poisson-Boltzmann theory [2j, to be further extended 
to the case of ubiquitous size-asymmetric dense ionic mix¬ 
tures. It is probably in this latter case that it will prove to 
be most useful specifically in the context of ionic liquids 

0 - 


For the Carnahan-Starling fluid we have found an 
asymptotic form that gives a rather simple analytic re¬ 
lation between potential and distance, eq. (51), as well 
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Figure 6: Top: calculation of the differential capacitance 
for a volume asymmetric lattice fluid. Bottom: volume 
asymmetric Carnahan-Starling fluid. From eq. (55). Each 
curve is for a different value of the chemical potential, 
p. Strongly negative chemical potentials give a minimum 
in the curves near ip = 0, together with two asymmetric 
maxima. As the chemical potential increases the curves 
develop a peak near (but not at) ip = 0. The chemical po¬ 
tential has been shifted such that the minimum of p{ip) is 
at ip = 0. In both cases the positively charged component 
has a volume 1.5 3 larger than the negative component. 


as the relation between potential and local charge density 
eq. (48). It is clear that the description of charged fluids as 
lattice gases or as charged hard spheres gives very differ¬ 
ent phenomenology in high field regions. The lattice gas 
crosses over very rapidly to a close packed system, whereas 
much higher fields are needed to compress the hard sphere 
system, leading to very slow cross-overs in 1 /ip 1 / 3 in phys¬ 
ical properties such as charge density. 


Appendices. — 

Numerical methods. We wrote numerical codes to 
study the double Legendre transformed free energy 


/(ci, c 2 ) - p(ci + c 2 ) - tH-ZiCi - z 2 c 2 ) (56) 


where ip = —ip. We do this by working with the effective 
coordinates 


n =(ci + c 2 ) 

q =(,zici - z 2 c 2 ) (57) 


So that we are interested in stationary points of the func¬ 
tion. 

f{n,q) - pn-ipq (58) 

where we have expressed the free energy as a function of 
the two independent coordinates, n the number density 
and q the charge density. 

We proceed by constructing an intermediate function 


Quid) by numerical minimisation of eq (58), with fixed p 
and q , with ip = 0. The function g^(q) is then passed 
to the Cliebfun library 22 23 which evaluates g^iq) for 


different specific values of q and builds a Chebyshev ap- 
proximant accurate to a relative accuracy of 10 _1 °. From 
this function we build the Legendre transform from q to 
ip by standard operations on g M 120 . 


9u (fl) -> flu (9) 


(flL) _1 W 


{g')-\iP')diP' 


(59) 


These steps are all performed by manipulation of the 
Chebyshev series, while maintaining close to machine pre¬ 
cision in the evaluations. The result is an approximant to 
p{ip). The last step is to transform back to p(ip) which 
requires a flip in sign of the potential axes. 

The functions p(ip) and g^{q) encode complementary in¬ 
formation on the physical system. We can find the equilib¬ 
rium charge density at a given potential from the relation 

= m 

we find the potential at imposed charge density from 

m = (si) 

aq 


The non-standard signs in these relations come from the 
difference between ip and ip. 

The question finally arrises as to how to use the numer¬ 
ically determined curves for p(ip) in other external codes. 
Inspiration comes from the Carnahan-Starling approxima¬ 
tion for the pressure which is a ratio of polynomials in the 
density. Such a general form is an example of a Pade ap¬ 
proximant that yields a high precision representation of 
the function p(ip) with an approximation as a ratio of two 
cubic polynomials that yields a rather good fit. Use of 
two quartics gives results which are visually perfect. Thus 
the present functional forms can be easily exported (this 
is even part of the chebfun library) to simple, fast approx¬ 
imations that can be used in other simulation codes. 

Clearly these methods are completely general can be 
applied to even more elaborate equations of state, extrap¬ 


olated from the best virial expansions 24 
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Convex formulation for Poisson-Boltzmann free en¬ 
ergies. As an alternative to writing the Poisson- 
Boltzmann functional in terms of the potential ij) with the 
help of the function p(ip) we can generate an equivalent 
convex formulation using the displacement field D. As 
shown in [25||26| this exact transformation requires the 
Legendre transfrom of the function p(if). However, we 
have already evaluated this object, it is just (q) , eq. (59). 
We can thus at once conclude that the general convex 
Poisson-Boltzmann function equivalent to those discussed 
above is 


D 2 


ff/i(div D - p e ) 


(62) 


with p e the external, imposed charged density. This form 
can be particularly interesting for the numerical work 
when coupling to other conformational degrees of freedom 
such as polymer chains or biomolecules. While we do not 
have analytic expression for g M for the Carnahan-Starling 
fluid it is again easy to generate the curve as a Chebyshev 
polynomial and export them to an accurate and efficient 
form for use in other codes. 
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